library(sp)
library(spdep)
library(dplyr)
library(rstan)
library(loo)
library(ggplot2)
library(bayesplot)
library(glmnet)

Error in library(glmnet) : there is no package called ‘glmnet’

GM@data <- GM@data %>% 
  mutate(
    log_pm25 = log(pm25), 
    log_sat = log(sat_2014)
  )

What variables are available in this data set?

print(colnames(GM@data))

[1] "City_locality"     "iso3"              "country"          
[4] "super_region"      "super_region_name" "pm25"             
[7] "sat_2014"          "log_pm25"          "log_sat"          

Regression on the original values instead of in log-domain

xylabs <- labs(
  x = expression(satellite), 
  y = expression(PM[2.5])
)
region_color_manual = scale_color_manual(
    values = 
      c("E-Eur/C-Eur/C-Asia" = "#00C094",
        "HighIncome" = "#FB61D7",
        "LatAm/Carib" = "#53B400",
        "N-Afr/MidEast" = "#A58AFF",
        "S-Asia" = "#00B6EB",
        "SE-Asia/E-Asia/Oceania" = "#C49A00",
        "Sub-Saharan Afr" = "#F8766D")
  )
plot1 <- ggplot(GM@data, aes(y = pm25, x = sat_2014)) +
  geom_point(aes(color = super_region_name), alpha = 0.4, size = rel(0.8)) +
  region_color_manual +
  geom_smooth(method = lm, color = "black", size = 0.5, linetype = 2) +
  coord_equal() +
  xylabs +
  guides(color = guide_legend(
    title = NULL,
    override.aes = list(alpha = 1, size = 2)
  )) +
  theme(legend.text = element_text(size = rel(0.6)))
plot(plot1)

Fitting a second degree polynomial instead

xylabs <- labs(
  x = expression(log(satellite)), 
  y = expression(log(PM[2.5]))
)
region_color_manual = scale_color_manual(
    values = 
      c("E-Eur/C-Eur/C-Asia" = "#00C094",
        "HighIncome" = "#FB61D7",
        "LatAm/Carib" = "#53B400",
        "N-Afr/MidEast" = "#A58AFF",
        "S-Asia" = "#00B6EB",
        "SE-Asia/E-Asia/Oceania" = "#C49A00",
        "Sub-Saharan Afr" = "#F8766D")
  )
plot1 <- ggplot(GM@data, aes(y = log_pm25, x = log_sat)) +
  geom_point(aes(color = super_region_name), alpha = 0.4, size = rel(0.8)) +
  region_color_manual +
  geom_smooth(method = glm, formula=y~poly(x, 2), color = "black", size = 0.5, linetype = 2) +
  coord_equal() +
  xylabs +
  guides(color = guide_legend(
    title = NULL,
    override.aes = list(alpha = 1, size = 2)
  )) +
  theme(legend.text = element_text(size = rel(0.6)))
plot(plot1)

Model used in the paper

xylabs <- labs(
  x = expression(log(satellite)), 
  y = expression(log(PM[2.5]))
)
region_color_manual = scale_color_manual(
    values = 
      c("E-Eur/C-Eur/C-Asia" = "#00C094",
        "HighIncome" = "#FB61D7",
        "LatAm/Carib" = "#53B400",
        "N-Afr/MidEast" = "#A58AFF",
        "S-Asia" = "#00B6EB",
        "SE-Asia/E-Asia/Oceania" = "#C49A00",
        "Sub-Saharan Afr" = "#F8766D")
  )
plot1 <- ggplot(GM@data, aes(y = log_pm25, x = log_sat)) +
  geom_point(aes(color = super_region_name), alpha = 0.4, size = rel(0.8)) +
  region_color_manual +
  geom_smooth(method = lm, color = "black", size = 0.5, linetype = 2) +
  coord_equal() +
  xylabs +
  guides(color = guide_legend(
    title = NULL,
    override.aes = list(alpha = 1, size = 2)
  )) +
  theme(legend.text = element_text(size = rel(0.6)))
plot(plot1)

r1 <- subset(GM@data, super_region==1)
r2 <- subset(GM@data, super_region==2)
r3 <- subset(GM@data, super_region==3)
r4 <- subset(GM@data, super_region==4)
r5 <- subset(GM@data, super_region==5)
r6 <- subset(GM@data, super_region==6)
r7 <- subset(GM@data, super_region==7)

b1 <- lm(r1$log_pm25 ~ r1$log_sat)
print(b1$coefficients)

(Intercept)  r1$log_sat 
  1.2928677   0.4628297 

b2 <- lm(r2$log_pm25 ~ r2$log_sat)
print(b2$coefficients)

(Intercept)  r2$log_sat 
   1.424214    0.700595 

b3 <- lm(r3$log_pm25 ~ r3$log_sat)
print(b3$coefficients)

(Intercept)  r3$log_sat 
  0.4585961   0.9131791 

b4 <- lm(r4$log_pm25 ~ r4$log_sat)
print(b4$coefficients)

(Intercept)  r4$log_sat 
  1.0401566   0.6807239 

b5 <- lm(r5$log_pm25 ~ r5$log_sat)
print(b5$coefficients)

(Intercept)  r5$log_sat 
  2.3130627   0.2960964 

b6 <- lm(r6$log_pm25 ~ r6$log_sat)
print(b6$coefficients)

(Intercept)  r6$log_sat 
  1.8787025   0.5071027 

b7 <- lm(r7$log_pm25 ~ r7$log_sat)
print(b7$coefficients)

(Intercept)  r7$log_sat 
  2.9589930   0.1482167 

plot_r1 <- ggplot(r3, aes(y = log_pm25, x = log_sat)) +
  geom_point(aes(colour = super_region_name), alpha = 0.2, size = rel(0.75)) + 
  region_color_manual +
  geom_smooth(aes(colour = super_region_name), method = lm) + 
  coord_equal() +
  xylabs + 
  legend_none() 
plot(plot_r1)

Same plot as above + separate fits for each super-region

plot2 <- ggplot(GM@data, aes(y = log_pm25, x = log_sat)) +
  geom_point(aes(colour = super_region_name), alpha = 0.2, size = rel(0.75)) + 
  geom_smooth(method = lm, color = "black", size = 0.5, linetype = 2) + 
  region_color_manual +
  geom_smooth(aes(colour = super_region_name), method = lm) + 
  coord_equal() +
  xylabs + 
  legend_none() 
plot(plot2)

Derive labels from the data based on average PM2.5 concentration in each country

average <- 
  GM@data %>% 
  group_by(iso3) %>% 
  summarise(pm25 = mean(pm25))
d <- dist(average)

NAs introduced by coercion

hh <- hclust(d)
clust <- cutree(hh,k = 6)
GM@data$cluster_region <-
  sapply(GM@data$iso3, function(x) clust[which(average$iso3 == x)])
print(table(GM@data$cluster_region))

   1    2    3    4    5    6 
  14 2172  400  376   10    8 

cluster_color_manual = scale_color_manual(
    values = 
      c("4" = "#00C094",
        "1" = "#FB61D7", 
        "5" = "#53B400",
        "2" = "#A58AFF",
        "3" = "#00B6EB",
        "6" = "#C49A00",
        "7" = "#F8766D")
  )

plot3 <-
  ggplot(GM@data, aes(
    y = log_pm25,
    x = log_sat
  )) + 
  geom_point(aes(colour = as.factor(cluster_region)), alpha = 0.2, size = rel(0.75)) + 
  cluster_color_manual +
  geom_smooth(
    method = lm, 
    color = "black", 
    size = 0.5, 
    linetype = 2
  ) + 
  geom_smooth(aes(colour = as.factor(cluster_region)), method = lm) + 
  coord_equal() +
  guides(color = guide_legend(
    title = NULL,
    override.aes = list(alpha = 1, size = 2)
  )) +
  xylabs

plot(plot3)

# Prior predictive simulations --------------------------------------------
# Inverse Gamme distribution. 
library(gridExtra)
library(grid)
library(invgamma)
library(ggplot2); theme_set(theme_bw())
x <- seq(0, 500, .01)
qplot(x, dinvgamma(x, 1, 1), geom = "line") -> plot1
qplot(x, dinvgamma(x, 1, 100), geom = "line") -> plot2
#qplot(x, 1/dinvgamma(x, 1, 1), geom = "line") -> plot3
#qplot(x, 1/dinvgamma(x, 1, 100), geom = "line") -> plot4
#qplot(x, dinvgamma(x, 2, 100), geom = "line") -> plot3
#qplot(x, dinvgamma(x, 4, 100), geom = "line") -> plot4
grid.arrange(plot1, plot2, ncol = 2)

# Plot: prior predictive with vague priors
set.seed(seed = 1923840483)
tau0 <- 1 / sqrt(rgamma(1, 1, rate = 100))
tau1 <- 1 / sqrt(rgamma(1, 1, rate = 100))
sigma <- 1 / sqrt(rgamma(1, 1, rate = 100))
beta0i <- rnorm(8, 0, tau0)
beta1i <- rnorm(8, 0, tau1)
beta0 <- rnorm(1, 0, 100)
beta1 <- rnorm(1, 0, 100)
Nsim <- length(GM@data$super_region)
xysim_labs <- labs(
  x = expression(paste("Observed ", log(PM[2.5]))),
  y = "Simulated data (log scale)"
)
data1 <- data.frame(
  log_pm25 = GM$log_pm25,
  sim = beta0 + beta0i[GM$super_region] +
    (beta1 + beta1i[GM$super_region]) * GM$log_sat +
    rnorm(Nsim, mean = 0, sd = sigma)
)
theme_set(bayesplot::theme_default(base_size = 18))
theme_update(axis.text = element_text(size = 20))
ggplot(data1, aes(x = log_pm25, y = sim)) + 
  geom_point(alpha = 0.1, color = "red") + 
  xysim_labs + ggtitle("Vague Recommended Priors")
ggsave(filename = "plots/prior_pred_vague.png", width = 4.5, height = 3.75)

# Plot: prior predictive with weakly informative priors
set.seed(seed = 1923840479)
tau0 <- abs(rnorm(1, 0, 1))
tau1 <- abs(rnorm(1, 0, 1))
sigma <- abs(rnorm(1, 0, 1))
beta0i <- rnorm(8, 0, tau0)
beta1i <- rnorm(8, 0, tau1)
beta0 <- rnorm(1, 0, 1)
beta1 <- rnorm(1, 1, 1)
data2 <- data.frame(
  log_pm25 = GM$log_pm25,
  sim = beta0 + beta0i[GM$super_region] +
    (beta1 + beta1i[GM$super_region]) * GM$log_sat +
    rnorm(Nsim, mean = 0, sd = sigma)
)
ggplot(data2, aes(x = log_pm25, y = sim)) +
  geom_point(alpha = 0.1) + 
  xysim_labs + ggtitle("Weakly informative joint prior data generating process")
ggsave(filename = "plots/prior_pred_wip.png", width = 4.5, height = 3.75)

# Plot: prior predictive with weakly informative priors (no seed)
tau0 <- abs(rnorm(1, 0, 1))
tau1 <- abs(rnorm(1, 0, 1))
sigma <- abs(rnorm(1, 0, 1))
beta0i <- rnorm(8, 0, tau0)
beta1i <- rnorm(8, 0, tau1)
beta0 <- rnorm(1, 0, 1)
beta1 <- rnorm(1, 1, 1)
data2 <- data.frame(
  log_pm25 = GM$log_pm25,
  sim = beta0 + beta0i[GM$super_region] +
    (beta1 + beta1i[GM$super_region]) * GM$log_sat +
    rnorm(Nsim, mean = 0, sd = sigma)
)
ggplot(data2, aes(x = log_pm25, y = sim)) +
  geom_point(alpha = 0.1) + 
  xysim_labs
ggsave(filename = "plots/prior_pred_wip.png", width = 4.5, height = 3.75)

# Plot: prior predictive comparison
data3 <- data.frame(
  log_pm25 = GM$log_pm25, 
  wip = data2$sim, 
  vague = data1$sim
)
ggplot(data3, aes(x=log_pm25, y=wip)) + 
  geom_point(alpha = 0.1) + 
  geom_point(
    aes(y = vague), 
    color = "red", 
    alpha = 0.1
  ) + 
  xysim_labs
ggsave(filename = "plots/prior_pred_compare.png", width = 4.5, height = 3.75)

# Fit Stan models 1, 2, 3 -------------------------------------------------
# Compile Stan programs
# * simple.stan: simple linear regression (Model 1)
# * hier.stan: non-centered parameterization of hierarchical model (Model 2, Model 3)
simple_mod <- stan_model("stan/simple.stan")
hier_mod <- stan_model("stan/hierarchical.stan")
# Data for model 1
standata1 <- with(GM@data, list(
  N = length(log_pm25),
  log_pm = log_pm25,
  log_sat = log_sat
))
# Data for model 2 (using super-regions from WHO)
standata2 <- with(GM@data, list(
  N = length(log_pm25),
  R = length(unique(super_region)),
  log_pm = log_pm25,
  log_sat = log_sat,
  region = super_region
))
# Data for model 3 (using super-regions from clustering)
standata3 <- standata2
standata3$R <- length(unique(GM$cluster_region))
standata3$region <- GM$cluster_region
# Fit the models with Stan
nuts_controls <- list(max_treedepth = 15, adapt_delta = 0.99)
mod1 <- sampling(simple_mod, data = standata1, seed = 2402)

SAMPLING FOR MODEL 'simple' NOW (CHAIN 1).

Gradient evaluation took 0.000554 seconds
1000 transitions using 10 leapfrog steps per transition would take 5.54 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)
Iteration:  600 / 2000 [ 30%]  (Warmup)
Iteration:  800 / 2000 [ 40%]  (Warmup)
Iteration: 1000 / 2000 [ 50%]  (Warmup)
Iteration: 1001 / 2000 [ 50%]  (Sampling)
Iteration: 1200 / 2000 [ 60%]  (Sampling)
Iteration: 1400 / 2000 [ 70%]  (Sampling)
Iteration: 1600 / 2000 [ 80%]  (Sampling)
Iteration: 1800 / 2000 [ 90%]  (Sampling)
Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 2.52046 seconds (Warm-up)
               2.69208 seconds (Sampling)
               5.21255 seconds (Total)


SAMPLING FOR MODEL 'simple' NOW (CHAIN 2).

Gradient evaluation took 0.000181 seconds
1000 transitions using 10 leapfrog steps per transition would take 1.81 seconds.
Adjust your expectations accordingly!


Iteration:    1 / 2000 [  0%]  (Warmup)
Iteration:  200 / 2000 [ 10%]  (Warmup)
Iteration:  400 / 2000 [ 20%]  (Warmup)

# Graphical posterior predictive checks -----------------------------------
theme_set(bayesplot::theme_default(base_size = 14))
y <- standata2$log_pm
yrep1 <- as.matrix(mod1, pars = "log_pm_rep")
yrep2 <- as.matrix(mod2, pars = "log_pm_rep")
yrep3 <- as.matrix(mod3, pars = "log_pm_rep")
samp100 <- sample(nrow(yrep1), 100)
# overlaid densities
color_scheme_set("blue")
ppc_dens_overlay(y, yrep1[samp100, ]) + 
  coord_cartesian(ylim = c(0, 0.7), xlim = c(0, 6)) +
  legend_none()
ggsave(filename = "plots/ppc_dens1.png", width = 4.5, height = 3.75)

color_scheme_set("gray")
ppc_dens_overlay(y, yrep2[samp100, ]) + 
  coord_cartesian(ylim = c(0, 0.7), xlim = c(0, 6)) +
  legend_none()
ggsave(filename = "plots/ppc_dens2.png", width = 4.5, height = 3.75)

color_scheme_set("red")
ppc_dens_overlay(y, yrep3[samp100, ]) + 
  coord_cartesian(ylim = c(0, 0.7), xlim = c(0, 6)) +
  legend_none()
ggsave(filename = "plots/ppc_dens3.png", width = 4.5, height = 3.75)

# stat: skew 
skew <- function(x) {
  xdev <- x - mean(x)
  n <- length(x)
  r <- sum(xdev^3) / sum(xdev^2)^1.5
  return(r * sqrt(n) * (1 - 1/n)^1.5)
}
color_scheme_set("blue")
ppc_stat(y, yrep1, stat = "skew", binwidth = 0.01) + 
  xlim(0, .6) + 
  legend_none()
ggsave(filename = "plots/ppc_skew1.png", width = 4.5, height = 3.75)

color_scheme_set("gray")
ppc_stat(y, yrep2, stat = "skew", binwidth = 0.01) + 
  xlim(0, .6) + 
  legend_none()
ggsave(filename = "plots/ppc_skew2.png", width = 4.5, height = 3.75)

color_scheme_set("red")
ppc_stat(y, yrep3, stat = "skew", binwidth = 0.01) + 
  xlim(0, .6) + 
  legend_none()
ggsave(filename = "plots/ppc_skew3.png", width = 4.5, height = 3.75)

# stat: group medians
superregion <- GM@data$super_region_name
superregion <- factor(superregion, levels = levels(superregion)[c(2,1,3:7)])
color_scheme_set("blue")
ppc_stat_grouped(y, yrep1, 
                 group = superregion, 
                 stat = "median", 
                 facet_args = list(nrow = 2)) + 
  facet_text(size = rel(0.7)) + 
  scale_x_continuous(breaks = function(x) pretty(x, n = 3)) +
  legend_none()
ggsave(filename = "plots/ppc_med_grouped1.png", height = 3, width = 7)

color_scheme_set("gray")
ppc_stat_grouped(y, yrep2, 
                 group = superregion, 
                 stat = "median",
                 facet_args = list(nrow = 2)) +
  facet_text(size = rel(0.7)) + 
  scale_x_continuous(breaks = function(x) pretty(x, n = 3)) +
  legend_none()
ggsave(filename = "plots/ppc_med_grouped2.png", height = 3, width = 7)

color_scheme_set("red")
ppc_stat_grouped(y, yrep3, 
                 group = standata3$region,
                 stat = "median", 
                 facet_args = list(nrow = 2)) + 
  facet_text(size = rel(0.7)) + 
  scale_x_continuous(breaks = function(x) pretty(x, n = 3)) +
  legend_none()
ggsave(filename = "plots/ppc_med_grouped3.png", height = 3, width = 7)